Numerical simulations of gravitational collapse in Einstein-aether theory 
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We study gravitational collapse of a spherically symmetric scalar field in Einstein-aether theory 
(general relativity coupled to a dynamical unit timelike vector field). The initial value formulation 
is developed, and numerical simulations are performed. The collapse produces regular, stationary 
black holes, as long as the aether coupling constants are not too large. For larger couplings a finite 
area singularity occurs. These results are shown to be consistent with the stationary solutions found 
previously. 
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I. INTRODUCTION 

Over the past several years there has been a renewal 
of interest in the physics of Lorentz symmetry violation. 
These investigations are largely motivated by the endur- 
ing problem of quantum gravity and also cosmological 
mysteries such as the dark matter and dark energy prob- 
lems. Lorentz violation (LV) effects in Standard Model 
extensions and modified particle dispersion relations are 
currently tightly constrained by high precision tests [H, 
however constraints on gravitational LV effects are not 
as severe. Various approaches to incorporating LV into 
gravity have been inve stig ated, see for example Refs. 

In this paper we will consider "Einstein-aether" theory 
(or ae-theory, for short) JEHU, which consists of a dy- 
namical unit timelike vector field u a coupled to Einstein 
gravity. The u a "aether" field can be thought of as the 
4-velocity of a preferred frame. The LV here is a form 
of spontaneous symmetry breaking, since the unit time- 
like constraint requires the vector field to have a non-zero 
value everywhere in any solution, even flat spacctimc. 

The Lagrangian of ae-theory depends on four dimen- 
sionless coupling constants 01,2,3,4- By studying the phe- 
nomenology of ae-theory it can be compared with obser- 
vations in order to constrain the c; values and possibly 
rule the theory out entirely as a description of nature. 
Observational constraints on ae-theory have been deter- 
mined from parameterized post-Newtonian analysis [HI 
[H HI, stability and energy_positivity [S H3, EE , pri- 
mordial nucleosynthesis [23| . and vacuum Cerenkov ra- 
diation [2l| . The combined constraints from all of these 
are reviewed in Ref. [Hj]. They are all consistent with a 
large (order unity) 2d region in the 4d parameter space. 
Agreement with weak field radiation damping (such as in 
binary pulsar systems) [3, HE} is expected to restrict the 



'Electronic address: garfinkl@oakland.edu 
^Electronic address: cteling@physics.umd.edu 
^Electronic address: jacobson@umd.edu 



parameters to a narrow band of width ~ 10~ 3 in that 2d 
region, and combined constraints from many such sys- 
tems may further narrow the allowed region to a small 
neighborhood of the origin, with size of order 10 -1 -10~ 2 . 

It is likely that theories with LV can generally be more 
tightly constrained by examining strong field solutions. 
In addition, the properties of such solutions may also be 
of interest from a purely mathematical physics perspec- 
tive since the additional fields that couple to GR can 
produce interesting new sorts of solutions and phenom- 
ena. Black hole solutions have been considered in the 
ghost condensate theory (26[, and spherically symmetric 
solutions in the relativistic MOND theory were studied 
in [27]] ■ Work on this subject in ae-theory was begun 
long ago in special cases jj, 0, @|, and a study of the 
solutions for general values of the coupling parameters 
was completed recently in [28| and (29j . 

In these papers the time independent spherically sym- 
metric solutions in ae-theory were characterized by ex- 
amining the field equations as a set of coupled ordinary 
differential equations. It was demonstrated that there 
is generically a three parameter family of solutions to 
these equations. Imposing the boundary condition of 
asymptotic flatness at spatial infinity reduces the solu- 
tion space to a two parameter family. This differs from 
GR, in which there is a unique spherically symmetric so- 
lution, the Schwarzschild metric, which is automatically 
asymptotically flat. The one parameter subset of solu- 
tions where u a is aligned with the timelike Killing vector 
was found analytically and discussed in detail in [281 ]. It 
was also shown that this solution describes the exterior 
of a static fluid star. 

Regular black holes in ae-theory must be described by 
a different solution, since the aether is a timelike vec- 
tor so it cannot coincide with the Killing vector which 
is null on the horizon and spacelike inside. That is, the 
aether must be flowing into the black hole. Moreover, 
the theory possesses spin-2, spin- 1, and spin-0 "massless" 
wave modes that travel at speeds that differ from one an- 
other and from the metric speed of light [T^] . To qualify 
as a black hole, a solution must therefore contain a re- 
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gion where all of these wave modes are trapped. It was 
found in [2t| that in spherical symmetry, regularity at 
the metric, spin-2, and spin-1 horizons is automatic, but 
regularity at the spin-0 horizon is a supplementary condi- 
tion. This condition reduces the two-parameter family of 
asymptotically flat solutions to a one-parameter family 
characterized by the total mass or horizon area. Particu- 
lar black hole solutions with regular spin-0 horizons were 
found by numerical integration for special cases of the 
coupling parameters. Like the Schwarzschild solution, 
these solutions have spacelike singularities at the origin, 
and in the exterior the metric functions are very close to 
the familiar Schwarzschild solution, even for fairly large 
coupling constants Cj. (The differences become more pro- 
nounced in the black hole interior.) For certain large val- 
ues of the coupling constants no black hole solution with 
a regular spin-0 horizon was found. 

From an astrophysical point of view, an important 
question is which solutions describe the final state of a 
gravitational collapse of matter fields. In [2^| it was con- 
jectured that the solutions with regular spin-0 horizons — 
when they exist — are the ones that arise from collapse of 
initially nonsingular data. The main objective of this 
paper is to test that conjecture by numerical evolution 
of the time dependent spherically symmetric field equa- 
tions given appropriate collapsing initial data. A mass- 
less scalar field coupled only to the metric will serve as 
the matter field. 

In Section fll] we review the nature and status of 
Einstein-aether theory, and discuss the nature of black 
holes in this theory. Section IIIII together with the Ap- 
pendix develops the initial value formulation under the 
assumption that the aether vector is hypersurface orthog- 
onal, as is always the case in spherical symmetry. The 
results of numerical collapse simulations are presented in 
Section ITVl We find that the collapsing scalar field does 
indeed reach a final stationary black hole state with a reg- 
ular spin-0 horizon, for the range of coupling constants 
found in [2^| to allow such states. Outside this range of 
Ci we find solutions that develop a finite area singularity 
at the would-be spin-0 horizon. The numerical evolution 
procedure also allows for wider explorations of space of 
Cj than was done in (29| . The results appear to support 
the conjecture that formation of black holes with regu- 
lar spin-0 horizons is generic in ae-theory for coupling 
parameters that are not too large. 



II. EINSTEIN-AETHER THEORY 

The action for Einstein-aether theory is the most gen- 
eral generally covariant functional of the spacctimc met- 
ric g a b and aether field u a involving no more than two 
derivatives (not including total derivatives), 



where 
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[R-K ab mn V a u m V b u n +X(g ab u a u b +l)} (2) 



and Lmatter denotes the matter lagrangian. Here R is the 
Ricci scalar, K ab mn is defined as 

K ab mn = c ig ab g mn + c 2 8 a m 8 b n + c 3 6^S b n ~ c iU a u b g mn (3) 

where the Ci are dimensionless coupling constants, and A 
is a Lagrange multiplier enforcing the unit timelike con- 
straint on the aether. The convention used in this paper 
for metric signature is ( — I — I — h) and the units are chosen 
so that the speed of light defined by the metric g a b is 
unity. In the weak-field, slow-motion limit ae-theory re- 
duces to Newtonian gravity [23[ , with a value of Newton's 
constant G*n related to the parameter G in the action (TIJ 
byGN = G , (l-( Cl +c 4 )/2)- 1 . 

In this paper we take the matter field to be a mini- 
mally coupled massless scalar field Xi with Lagrangian 
— V a xV a x- However, to simplify the form of the equa- 
tions of motion, we introduce the quantity ip = \ vl67rG 
so that the matter Lagrangian takes the form 
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(4) 



S = 



(1) 



The field equations from varying (fT]) with respect to 
g ab , u a , ip, and A are given respectively by 

Gab — T a b (5) 

V a J a b + Xu b + c 4 a Q V h u a = (6) 
V a V a ^ = (7) 

U a U a = -1. (8) 

Here G a b is the Einstein tensor of the metric g a b] and the 
quantities J a b, a a and the stress-energy T a b are given by 

J a m = K ab mn V b u n (9) 
a a = u b V b u a (10) 

Tab = -yab(j C dV c U d + V c ^V C ij) 

+ V a ^V 6 ^ + c 4 a a a b + Xu a u b 
+ a {V a u c V b u c - V c it a V c ii 6 ) 

+ V c [j C ( a M 6) + U C J [ab ) - J(a C U b )] ■ (H) 

In this paper we will perform numerical simulations 
of this system in spherical symmetry. The condition 
of spherical symmetry implies that the twist uj a — 
£abcdU b V c u d of the vector field u a vanishes. This makes 
the four parameters of the Lagrangian redundant: the 
action depends only on C2, C13 and Ci 4 where we use the 
abbreviation C13 to stand for c\ + C3 (and correspond- 
ingly for C14). (In Ref. [29[ it is shown that when the 
twist vanishes the C4 term can absorbed by the shifts 
ci — > ci + C4 and C3 — > C3 — C4. Thus the action depends 
only on C2, C14 and C3 — C4. But since the latter is equal 
to C13 — C14, the action is also determined by C2, C13, and 
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ci4.) Because the twist vanishes, it follows that u a is 
hypersurface orthogonal. 

For our study of scalar held collapse we considered a 
case where the parameters Cj satisfy the requirements 
that (i) the Einstein-aether theory has precisely the same 
post-Newtonian parameters as general relativity, so it 
agrees in the post-Newtonian approximation with all ex- 
perimental tests of general relativity, (ii) the theory is 
linearly stable, (iii) the wave modes all carry positive en- 
ergy, and (iv) the wave modes travel at greater than or 
equal to the matter speed of light, so there is no vacuum 
Cerenkov radiation. Together these impose the condi- 
tions QJ] 



C>2 



2ci + c 3 



3ci 



ci 

< C13 < 1 

< (ci - c 3 ) < 



Cl3 



3(1 -da) 



(12) 

(13) 
(14) 
(15) 



We chose c\ = 1/3 and C3 = 1/6 (essentially in the 
midrange of the inequalities given above) and then de- 
termined C2 and C4 by the above equalities. This yields 
the parameter set 



1 

Cl= 3' 



C2 



1 

T C3 



-, c 4 = — -. (16) 
6' 12 y ' 



This choice of the constants Cj is different from those 
used in the studies of static black holes in [29(. However, 
for the Einstein-aether theory without matter, the trans- 
formation u a — * \fou a (with a a constant and with h a b 
remaining unchanged) takes Einstein-aether theory into 
itself but with changes in the parameters c; [30]. This 
transformation and the restriction to spherical symmetry 
allows one to transform the results with our parameter 
choice (fT6"l) to one with C13 and C4 vanishing and with 
ci = 1/4 and C2 = 1/2. This new set of parameters falls 
in the class that lends itself to a numerical integration of 
the static, spherically symmetric Einstein-aether equa- 
tions using the methods of [29(, and thus allows a direct 
comparison of a static solution with the end state of our 
numerical simulations of gravitational collapse. 

In this paper we examine the collapse process to de- 
termine the nature of the black hole that is formed. 
This issue is somewhat more subtle in Einstein-aether 
theory than in general relativity. Usually in a gravita- 
tional collapse simulation one looks for apparent hori- 
zons: marginally outer trapped surfaces of the metric 
g a b- Typically the spacetime outside the apparent hori- 
zon settles down to a stationary (static in the spherically 
symmetric case) black hole. However, Einstein-aether 
theory contains a spin-0 mode that travels at the speed 
vq where fijl 



ci 23 (2 - c i4 ) 



ci 4 (l - ci 3 )(2 + ci 3 + 3c 2 ) ' 



(17) 



For the equalities that we have imposed on the Ci this 
expression reduces to 
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3(d - c 3 )(l - ci 3 ) ' 



(18) 



Thus, the boundary of the black hole should be a horizon 
that just barely traps the spin-0 mode. This turns out 
to be equivalent to a marginally outer trapped surface of 
the metric g a b = <? a b + (l — v$)u a Ub- Thus in searching for 
black holes we look for marginally outer trapped surfaces 
of the metric g a b, which for brevity we will call spin- 
horizons. For the parameter choice (|TB)) . the spin-0 
speed from (fl~8|) is \/2, so the spin-0 horizon lies inside 
the metric horizon. 



III. INITIAL VALUE FORMULATION 

Numerical evolution requires a choice of the surfaces of 
constant time, and we will choose those surfaces orthog- 
onal to u a . This choice is always possible in spherical 
symmetry since u a is necessarily hypersurface orthogo- 
nal. The spatial metric h a b and extrinsic curvature K a b 
of these hypersurfaces arc then given by 



h a b - 
K ab 



g a b + u a u b 

_ _lr u 
~ n*^u ll ab 



(19) 
(20) 



where C denotes Lie derivative. Notice that the con- 
straint §8§ that u a be a unit vector is built into our ini- 
tial value formulation by its identification with the unit 
normal. 

The use of K a b as a variable allows us to write Ein- 
stein's equation in a form that is first order in time. Sim- 
ilarly we introduce the quantity P = C u ip which allows 
us to do the same for the wave equation. 

The time evolution vector held takes the form t a = 
au a + (3 a . We use as a radial coordinate r the length 
in the radial direction (rather than an area coordinate). 
This makes the spatial metric 



h ab = d a rd b r + $ 2 Y, a b 



(21) 



where S a b is the unit two-sphere metric. 

The quantities that are evolved are (ip, P, K, a r , $). All 
other quantities are found from these. In the appendix, 
equations of motion for these quantities are derived from 
the Euler-Lagrange equations (0[8]) . These equations are 



d t tp = aP + p r d r ip 
d t P = (3 r d r P 



a 



PK + a r d r ip + drdrijj + —^—drtp 



(22) 
(23) 



d t K = (3 r d r K- 
a 



a 



K 



2 + d3 + 3c 2 



(d4 — 2)(d r a r + 2a r 9 r $/$ + a r a r ) 
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2P 2 + 3(1 -c 13 )Q 2 



dta r = (3 r d r a r - 

Cl3 



-Pd r lp 



C123 



dt<S> = 



cu(l-c 13 ) c 14 (l-c 13 ) 
(3 r d r <P + a<f>(Q/2- K/3). 



d r K 



(24) 

(25) 
(26) 



Here Q = K r r — (K/3) is the trace-free part of the extrin- 
sic curvature. As shown in the appendix, the quantities 
Q, a and (3 are determined by 



d r Q 



3Q 
3> 



9,.$ + (1 - dsK 1 [i(2 + c 13 + 3c 2 )d r K 



d r In a — a r 
8 r f3 r = a(Q + K/3). 



(27) 
(28) 
(29) 



The Hamiltonian initial value constraint is 

>2 



C = d r d r $ 



1 



2$ 



Ci4a r 9,<i> 



cu(2d r a r + a r a r ) + P 2 + (d r ipy 



+ |(l-c 13 )g 2 -i(2 + c 1 3 + 3c 2 )if 5 



0. 



(30) 



For initial data, we choose a moment of time symmetry, 
so that P, Q and K are zero. Both a r and ip can be freely 
specified. We choose a r to vanish and ip to be 



ip = a,Q exp [- 



(31) 



where ao, tq and s are constants. Thus the scalar field 
is initially at rest and is essentially a spherical shell of 
radius ro, thickness s and amplitude ao. Then equa- 
tion (|3H|) is integrated outwards to find <&, subject to the 
boundary conditions $ = and <9 r $ = 1 at r = 0. The 
other boundary conditions that need to be imposed are 
(3 r = and Q = for smoothness at r = 0, and a — > 1 
for a gauge fixing as r — > oo. 

The evolution proceeds as follows: given 
(i/j, P, K,a r ,&) at one time step, equations (|28|) 
and (|27[) are integrated to find a and <5- Then equation 
([2"9"]1 is integrated to find /3 r . Then equations ([2"2"]l. 
(|2"3"|) . (J2U), (123 and (J25]) are used to find respectively 
(ip, P, K, a r , 4>) at the next time step. Since $ is evolved 
using equation (|26p but must also satisfy equation (|30|) 
this means that equation (|3H)l can be used as a code 
check in a convergence test. 



IV. NUMERICAL RESULTS 

The specific parameters we used for the initial scalar 
field are ao = 0.15, tq = 10 and s = 4. Initial data with 
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FIG. 1: Plot of C vs r for a simulation with 12800 points (solid 
line) and C/4 for a simulation with 6400 points (dashed line). 
The agreement of the two curves demonstrates second order 
convergence of the code. The curves are for a time when the 
scalar field is in the process of collapsing. 



these parameters leads to gravitational collapse and the 
formation of a black hole. 

In the numerical simulations spatial derivatives were 
implemented with standard centered finite differences. 
The time evolution was done with the Iterated Crank- 
Nicholson (ICN) method. Kreiss-Oliger type dissipation 
was used for stability. The first grid point was at r = 0. 
Of the quantities that are evolved by the code, smooth- 
ness requires that a r and $ are odd functions of r and 
that the other quantities are even. The code evolves the 
quanties on all grid points except the first and last. The 
odd quantities are set to zero on the first grid point, while 
the even quantities are set as follows: 



fx = (4/ 2 - / 3 )/3 



(32) 



which implements the condition that the derivative of / 
vanish at the origin. Here / denotes any even quantity, 
and the subscript denotes the corresponding gridpoint. 
At the last gridpoint, N, all the evolved quantities are 
set to zero with the exception of $ which is set by 



N 



2$ W _i - $tv-2 



(33) 



which is the appropriate boundary condition for a quan- 
tity that asymptotically satisfies d^/dr = 1. 

The simulations were run on a SunBlade 2000 in double 
precision. The number N of gridpoints used was 12800 
and the initial value of r at the outermost gridpoint was 
80. Figure |T]) shows the result of a convergence test 
of the code. It is a plot of the constraint quantity C 
of equation (|30|) done at a time when the scalar field is 
in the process of collapsing. The solid curve represents 
C for a simulation with 12800 points, while the dashed 
curve represents C/4 for a simulation with 6400 points. 
The fact that the two curves argee demonstrates that the 
code is second order convergent. 

The simulations show that for the initial data used the 
collapse does form a spin-0 horizon (which is contained 
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FIG. 2: Plot of a r vs z for the end state of the collapse sim- 
ulation (solid line) and the static solution found using the 
method of [2tJ (dashed line) 
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FIG. 4: Plot of Q v s z for the end state of the collapse sim- 
ulation (solid line) and the static solution found using the 
method of [2tJ (dashed line) 





FIG. 3: Plot of K vs z for the end state of the collapse sim- 
ulation (solid line) and the static solution found using the 
method of [29l ] (dashed line) 



FIG. 5: Plot of V a $V a $ vs z for the end state of the collapse 
simulation (solid line) and the static solution found using the 
method of [29[ (dashed line) 



within an ordinary horizon, that is a marginally outer 
trapped surface of the metric g a b)- Furthermore, outside 
of the spin-0 horizon the metric and aether field do settle 
down into a time independent state. In addition, all the 
scalar field either falls into the black hole or escapes, 
so that near the horizon the time independent state is 
that of a pure Einstein-aether black hole with no matter 
other than the aether field u a . We are able to follow the 
time evolution to late times with our horizon crossing 
time slices, orthogonal to the timelike aether field u a , 
because the lapse a is driven to zero as the singularity is 
approached. The reason for this can be traced to Eqn. 
(|28f . Evidently at the singularity the aether acceleration 
component a r goes to positive infinity. 

In reference 29] the static black holes of Einstein- 
aether theory with a regular spin-0 horizon were found 
numerically, and it was conjectured that these would be 
the endstates of gravitational collapse. Our simulations 
verify that this is indeed what happens. In figures (|2j). 
([3]) and Q are plotted respectively the quantities a r , K 
and Q for this simulation (solid line) and for the static so- 



lution found using the method of [29| (dashed line). Here 
the static solution is found using the transformed choice 
of the parameters c^, and its metric and aether variables 
are then transformed back as explained in Section [TT] to 
make the comparison with the result of the collapse simu- 
lation. In all cases the quantities are plotted as functions 
of z = $ — <&o where $o is the value of $ at the spin-0 
horizon. Note that the agreement is extremely good, so 
the end result of collapse is indeed the static black hole 
solution found using the method of [29| . 

The quantities a r , K and Q completely determine the 
aether field, and since we are treating spherically sym- 
metric spacetimes where there are no separate gravi- 
tational degrees of freedom, these quantities also com- 
pletely determine the metric. Nonetheless, it is helpful 
to do a direct comparison of the metric of the endstate of 
the collapse simulation and the static metric found using 
the method of 29]. For a static, spherically symmetric 
metric with area radius <I>, the metric is determined by 
the quantities V Q <I>V a <I> and V a V a <I>. In figures © and 
© these quantities are respectively plotted. As with 
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FIG. 6: Plot of V a V a $ vs z for the end state of the collapse FIG. 7: Plot of K vs r for simulations with ci = 0.7 (solid 
simulation (solid line) and the static solution found using the line) and ci = 0.8 (dashed line), 
method of [29j (dashed line) 



figures (J2EJ) the quantities are plotted as functions of z 
with the endstate of the collapse simulation as a solid line 
and the static metric found using the method of reference 
[291 ] as a dashed line. The agreement between the curves 
shows that the metric in the final state of the collapse 
simulation agrees with the static metric. 

We now present the results of a slightly different sim- 
ulation that uses directly a different parameter choice of 
[29l | . namely C3 = C4 = 0, and c-i chosen so that the spin-0 
horizon coincides with the usual metric horizon. That is, 
the speed vq in (jTTJ) is set equal to one, which requires 



(>2 



2 - 4ci + 3cf 



(34) 



If < c\ < 1 such theories meet all the linear stability 
conditions (stable modes, positive energy, and no vacuum 
Cerenkov radiation). In this class the post-Newtonian 
preferred frame parameter a\ vanishes identically, but ai 
does not vanish unless c\ — 2/3 (see [13, EH for 0:1,2)- 
Nevertheless, since the stability conditions are met, it 
seems likely that the behavior of collapse in the theory is 
qualitatively similar to the previous case where also ai 
vanishes. In [29(, static, asymptotically flat black hole so- 
lutions were found in the range c± < 0.7, by integrating 
from a regular horizon to spatial infinity using a shoot- 
ing method. The asymptotically flat boundary condition 
could not be met for c\ > 0.8. 

We performed collapse simulations for the case c\ = 0.7 
and the case c\ = 0.8. In the c\ — 0.7 case a regular 
horizon forms. However, in the c\ = 0.8 case, no horizon 
forms and the evolution seems to become singular, thus 
indicating the formation of a naked singularity. Figure 
([7]) Shows the value of K (dashed line) for the ci = 0.8 
simulation at a time shortly before the singularity forms. 
For comparison, the same figure shows (solid line) the 
value of K at the same time for the c\ = 0.7 simulation. 
Since the collapse solutions are by construction asymp- 
totically flat, this is consistent with what was found for 
static solutions in 12911. 



V. CONCLUDING REMARKS 

We have examined the process of dynamical collapse of 
a spherically symmetric, self-gravitating scalar field pulse 
in Einstein-aether theory. In pure Einstein gravity it is 
well known that for appropriate initial conditions such a 
pulse can form a black hole. Ours is the first study of 
the analogous process in Einstein-aether theory. In or- 
der to carry out the numerical simulations we developed 
the initial value formulation for this theory in spherical 
symmetry. 

Here we have minimized the initial "distortion" of the 
aether field, by setting its divergence and acceleration to 
zero. (At the other extreme, one could study collapse 
of a pure spherical aether pulse, with no matter field at 
all.) These initial conditions are not precisely what would 
occur in an astrophysical setting, but should be indica- 
tive. We found that for the initial data studied, black 
holes with regular metric and spin-0 mode horizons form 
when the parameters in the Lagrangian are not too large. 
For larger parameters, a singularity occurs at finite area 
where the spin-0 horizon would have been. These results 
are consistent with findings from previous investigation 
of static solutions [2!| . 

Since regular black holes have almost certainly been 
"observed" in nature, the black hole formation found 
here indicates that a strong-field test has been passed 
by ae-theory, at least for coupling parameters that are 
not so large that the spin-0 horizon is singular. Con- 
straints from the absence of vacuum Cerenkov radiation 
require that the spin-0 mode travels at a speed greater 
than or equal to the "speed of light" as seen by matter 
fields [2l[, here the "metric" speed of light. Thus a the 
singular spin-0 horizon would necessarily be on or inside 
the matter horizon. Whether it could then have any ob- 
servational consequences remains to be determined. 

It would also be interesting and astrophysically im- 
portant to examine the case of collapse with angular 
momentum, to a spinning black hole. The stationary 
spinning solutions have not yet been studied however, 
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and to follow the time-dependent collapse would require 
the general initial value formulation, without assuming 
hypersurface orthogonal aether nor spherical symmetry. 
In addition, axisymmetric codes are much more compli- 
cated. 
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APPENDIX: DERIVATION OF THE INITIAL 
VALUE FORMULATION FOR SPHERICALLY 
SYMMETRIC EINSTEIN-AETHER THEORY 

In this section we derive the evolution equations for 
Einstein-aether theory in spherical symmetry. To begin 
with, we use only the hypersurface orthogonality of u a 
and the Euler-Lagrange equations of motion ([5][5]). 

Equation ([6]) (the equation of motion for the aether 
field) then becomes 

cm (£ u a a + 2K ab a b - Ka a ) + c 13 D b K ab + c 2 D a K = 0. 

(A.l) 

Here C denotes the Lie derivative and D a is the spatial 
covariant derivative. 

The usual initial value formulation for general relativ- 
ity requires that we evolve h ab and K ab as well as the 
matter fields. From equation (|20|) we have 



(A.2) 



which is used to evolve h ab . The quantity K (the trace 
of K a b ) sasfies the evolution equation 

C U K = -D a a a - a a a a + K ab K ab + \T ab (h ab + u a u b ). 

(A.3) 

In addition, the quantity K ab satisfies the constraint 
equation 

D b K ab - D a K = -h b u c T bc . (A.4) 

The extrinsic curvature K a b can be decomposed into a 
trace and a trace-free part. We will use equation (|A.3[) 
to evolve the trace K and will use equation (|A.4[) to find 
the trace- free part. 

However, there is a subtlety associated with the 
Einstein-aether theory that makes the initial value prob- 
lem somewhat complicated: the stress energy tensor con- 
tains terms that have second time derivatives of the 
aether field and second time derivatives of the metric. 
Thus after writing out the Einstein field equations one 
must solve them for these second time derivatives. In 



particular, some straightforward but tedious algebra ap- 
plied to equation yields 

T ab (h ab + u a u b ) = -(c 13 + 3c 2 )£ u K + c u (D a a a + a a a a ) 
+ 2P 2 - 2c 13 K ab K ab + c 123 K 2 (A.5) 
-PD a iP 

d4 (2K ab a b - Ka a + C u a a ) (A.6) 



-h b u c T bc 



where P — C u ^ and c\ 23 — c\ + c 2 + c 3 . Using equation 
(|A.5[) in equation (|A.3|) and solving for C U K yields 

(2 + c 13 + 3c 2 ) C U K = (c 14 - 2){D a a a + a a a a ) + 2P 2 
+ 2(1 - c 13 )K ab K ab + c 123 K 2 . (A.7) 



Similarly, equation (IA.4|) . upon substituting the expres- 
sion in equation (|A.6[) , does not directly yield a constraint 
since the right hand side of equation (|A.6|) contains a 
term proportional to C u a a - However, eliminating this 
term using equation (| A. 1 [) yields the following constraint 

[1 - c 13 ]D b K ab = (1 + c 2 )D a K - PD a ^j (A.8) 

which can be used to rewrite the equation of motion for 

a a as 

jO u a a = -2K ab a b + Ka a H 77-^ -P£) Q ^ 



cu{l ~ C13) 



C123 



-D a K. 



(A.9) 



cu(i - ci 3 y 

The usual Hamiltonian constraint of general relativity 



(3) i? + K 2 - K ab K ab = 2T ab u a i 



(A.10) 



where is the scalar curvature of the spatial metric. 
However, it follows from equation (jTTJ) that 

2T ab u a u b = 2c m D a a a + P 2 + D a iJjD a iP 

+ c 14 a a a a - c 2 K 2 - c 13 K ab K ab . (A.ll) 



On substituting equation (|A.11|) into equation (|A.10[) one 
finds the constraint 

= c 14 (2D a a a + a a a a )+P 2 + D a 4;D a 4; 

+ (l-c 13 )K ab K ab ~(l + c 2 )K 2 . (A.12) 

We now impose spherical symmetry. We use as a radial 
coordinate r the length in the radial direction (rather 
than an area coordinate). This makes the spatial metric 

Kb - CVCV + $ 2 Sa6 (A.13) 

where T, ab is the unit two-sphere metric. The time evo- 
lution vector field takes the form t a = au a + j3 a . From 
the definition of P we find 



(A.14) 
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while from the wave equation we obtain 
dtP = (3 r d r P + a 



PK + a r d r ij) + d r d r ij) + ' d r ip 



(A.15) 

It is helpful to define Q = K r r — (K/3) so that Q is the 
component of the trace- free part of K a f,. From equation 
(IA.2[) , we have the standard result 



C-thab = —2aK ab + Cphab. 
The rr component of equation (|A.16j) yields 

d r (3 r = a(Q + K/3) 
while the 99 component yields 

= /T<9 r $ + a$(Q/2 - K/3). 
From the definition of a we have 
d r In a = a r ■ 

Equation (|A.8|) yields 



(A.16) 
(A.17) 
(A.18) 
(A.19) 



3Q 



while equations (|A.9|) and (| A.7|) yield respectively 



dta r = (3 r d r a r + a 

Cl3 



( 2K n\ 
{— ~ Q)a r 



P8 r 1p 



C123 



Ci 4 (l - C13) 

d t K = [3 r d r K+^K 2 



Cu(l - C13) 



d r K 



(A.21) 



2 + ci 3 + 3c 2 
2P 2 + 3(1 -c 13 )Q 2 



(C14 — 2)(d r a r + 2a r 9 r $/$ + a r a r ) 

(A.22) 



The Hamiltonian constraint (equation (|A.12[) ) becomes 
the vanishing of the quantity C where 



C = d r d r ® 



2$ 



Ci4a r d r $ 



cu(2d r a r + a r a r ) + P 2 + (d r ipY 



d r Q = -^9 r $+(l - ci 3 ) 1 [|(2 + C13 + 3c 2 )d r K - Pd r ^] 



+ |(l-c 13 )Q 2 -i(2 + c 1 3 + 3c 2 )if 2 



. (A.23) 



(A.20) 
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